module Qmat_Mod          
  use CommData
  use Phi_Mod
  use Operations_Mod
  implicit none
contains
!************************************************************************************
    subroutine Qm(k,i,Matr_Q)
        real*8,Intent(in)  :: k(:)
	    real*8,Intent(out) :: Matr_Q(:,:)
	    integer i,j
	    real*8 x,y,z,a,b,c
	    Matr_Q = 0.0
        Matr_Q(1,5) = Phi_kk(1,2,3,k,i)
        Matr_Q(1,9) = Phi_kk(1,3,2,k,i)
        Matr_Q(5,9) = Phi_kk(2,3,1,k,i)
	    x = Phi_kmk(1,2,3,k,i)
	    y = Phi_kmk(3,1,2,k,i)
	    z = Phi_kmk(2,3,1,k,i)
	    a = Phi_kpk(1,2,3,k,i)
	    b = Phi_kpk(3,1,2,k,i)
	    c = Phi_kpk(2,3,1,k,i)
        Matr_Q(2,4) = 0.5d0 * (x - a)
        Matr_Q(3,7) = 0.5d0 * (y - b)
        Matr_Q(6,8) = 0.5d0 * (z - c)
        Matr_Q = Matr_Q + transpose(Matr_Q)
        Matr_Q(1,1) = Phi_kk(1,1,3,k,i)
        Matr_Q(2,2) = 0.5d0*(x+a)
        Matr_Q(3,3) = 0.5d0*(y+b)
        Matr_Q(4,4) = 0.5d0*(x+a)
        Matr_Q(5,5) = Phi_kk(2,2,3,k,i)
        Matr_Q(6,6) = 0.5d0*(z+c)
        Matr_Q(7,7) = 0.5d0*(y+b)
        Matr_Q(8,8) = 0.5d0*(z+c)
        Matr_Q(9,9) = Phi_kk(3,3,1,k,i)
    end subroutine Qm
    !************************************************************************************
    subroutine Qmatrix(Eff, Plast, Qreal, i)
        real*8,intent(inout) :: Qreal(:,:)
	    real*8,intent(in) :: Eff(:),Plast(:)
	    INTEGER, intent(in) :: i
	    real*8  UU(9,9), VV(9,9), QQ(9,9), BB(9,9), II(3,3)
        real*8 A(3,3),U(3,3),VT(3,3),sing(3)
        real*8 work(20)
        integer lwork,info 
        lwork = 20
        II = 0d0
        II(1,1) = 1d0; II(2,2) = 1d0;II(3,3) = 1d0;
        A=reshape(Eff,(/3,3/))
        UU=0d0
        VV=0d0   
        call dgesvd('A','A',3,3,A,3,sing,U,3,VT,3,work,lwork,info)				
	    UU(1:3,1:3)=U;
	    UU(4:6,4:6)=U;
	    UU(7:9,7:9)=U;
	    call KronProd(transpose(VT), II, VV)
	    call Qm(sing,i,QQ)
	    call KronProd(reshape(Plast(:),(/3,3/)),II,BB)
	    QQ=matmul(BB,matmul(matmul(UU,matmul(VV,matmul(QQ,matmul(transpose(VV),transpose(UU))))),transpose(BB)))
	    Qreal(:,:) = QQ
    end subroutine Qmatrix
    !************************************************************************************
    subroutine Qmatrix_new(QQ, rotate_matrix, Qx)
        real*8,intent(in) :: rotate_matrix(:,:)
        real*8,intent(out) :: QQ(:,:), Qx(:,:)
        real*8 :: R(6,6), I(3,3)
        I = 0d0; R = 0d0;
        I(1,1) = 1d0; I(2,2) = 1d0;I(3,3) = 1d0;
        call KronProd(rotate_matrix,I,R)
        QQ = matmul(matmul(transpose(R),QQ(1:6,1:6)),R)
        Qx(:,:)= (QQ(1:3,1:3)+transpose(QQ(1:3,1:3)))/2d0
    end subroutine Qmatrix_new
end module Qmat_Mod